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Abstract 

We investigate the properties of strongly interacting bosons in two dimensions 
at zero temperature using mean-field theory, a variational Ansatz for the 
ground state wave function, and Monte Carlo methods. With on-site and 
short-range interactions a rich phase diagram is obtained. Apart from the 
homogeneous superfluid and Mott-insulating phases, inhomogeneous charge- 
density wave phases appear, that are stabilized by the finite-range interaction. 
Furthermore, our analysis demonstrates the existence of a supersolid phase, 
in which both long-range order (related to the charge-density wave) and off- 
diagonal long-range order coexist. We also obtain the critical exponents for 

the various phase transitions. 
PACS numbers: 67.90,+z, 05.30.Jp 
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I. INTRODUCTION 



In recent years strongly interacting bosons attracted a lot of attention. Experiments on 
thin granular films show a zero temperature superconductor to insulator (SI) phase transition 
as a function of the film thickness [|I| or magnetic field 0, and experiments on Josephson 
junction arrays (JJA) exhibit the same SI transition as a function of the ratio between 
Josephson coupling and charging energy ||. These findings are interpreted in terms of a 
bosonic model in which the competition between the hopping of Cooper pairs and the 
Coulomb interaction between them is responsible for the SI transition. In relation to the 
effect of disorder on the SI transition a Bose-glass phase is discussed in Refs. 

An interesting extension of the minimal model with on-site interaction only is to account 
for finite-range interactions, which are present in Josephson junction arrays and 4 He-films 
on substrates in two dimensions, and bulk 4 He in three dimensions. Artificially fabricated 
JJA's are ideal model systems to study the superconductor to Mott-insulator transition as a 
function of the coupling constants or magnetic field ||. The finite range of the interaction in 
these systems leads to commensurability and frustration effects. Furthermore, coexistence 
phases may appear that are called supersolids. In these phases a charge-density wave (CDW) 
that is stabilized by the interaction may coexist with superfluidity. In other words the system 
has diagonal long-range order (LRO) and off-diagonal LRO at the same time. 

The supersolid phase was studied in the early seventies after Andreev and Lifshitz || 
suggested that vacancies in a quantum crystal such as 4 He might Bose-Einstein condense 
at low temperatures without destroying the crystal structure, thereby establishing a super- 



fluid solid JT0[, or supersolid. Normally bosons at zero temperature are either superfluid 
(with off-diagonal LRO) or solid (with diagonal LRO). However, for a finite range of the 
interactions between the bosons a coexistence phase was predicted within mean-field approx- 
imations [11-14]]. Experiments have been performed on 4 He, but no positive identification of 
this coexistence phase has yet been made. There are, however, hints towards such a phase 

MM- 



2 



Recently several other kinds of coexistence phases were studied. The possibility of a 
spontaneous vortex anti-vortex lattice in superfluid films was explored in Ref. [|I7J] and a 



coexistence phase of superfluidity and hexatic orientational order was proposed in Ref. [18 



Orientational order in incompressible quantum Hall fluids is discussed in Ref. ||19|| . Collinear 
supersolids were studied in Refs. PU| , [2~T| . Finally, we mention the relation between 2D bosons 
and 3D flux-lines in type II superconductors (high T c materials) in a magnetic field |22H2"4|1 ■ 
Also in these systems different kinds of LRO may coexist and the equivalent of the supersolid 



is discussed in Refs. p5|j26|] . Related is the question whether or not vortices may form a 
disentangled liquid, which would imply a normal ground state for bosons with long-range 



Coulomb interaction [24 



In this paper both the SI transition and the supersolid phase will be studied in some 
detail. In the next section several related models for interacting bosons will be discussed, of 
which the Bose-Hubbard model is the most general. In the limit of infinite on-site interaction 
the Bose-Hubbard model reduces to the spin | XXZ model. For a large number of bosons 
per lattice site the Quantum-Phase model is applicable. The phase diagram of the latter two 
models are first determined within mean-field theory in Section |ITT[ A variational Ansatz for 
the ground state wave function that treats all three models on equal footing [27] is presented 
in Section |iy]. In Section [V| more accurate results for the phase diagram will be obtained by 
means of Quantum Monte Carlo. Part of the data has already been published as Ref. p8| . 
Finally the paper closes with a discussion of the various results in Section |VI|. 



II. THREE MODEL HAMILTON! ANS 

A convenient starting point for the description of interacting lattice bosons is the Bose- 
Hubbard Hamiltonian 

Hbh = \Y, n i u ii n i ~ V n i ~ | Y,( b l b i + h ) h i) > i 1 ) 
z ij i A (ij) 

where b\b are the creation and annihilation operators for bosons that satisfy the commu- 
tation relation [bi, fej] = 5ij and rij = b\bi is the number operator. The first term describes 
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the density-density interaction between the bosons. We will take it to be short range, Uq 
for on-site, U\ for nearest neighbors, and U-z for second-nearest neighbors. The second term 
describes the coupling of the density to a chemical potential, and the third describes the 
hopping from site to site with hopping integral t. The first two terms will cause the bosons 
to be localized if the interaction dominates over the hopping. In the opposite case where 
the hopping dominates the bosons will form a superfluid. 

In the limit of large on-site interaction Uq, only the states with and 1 boson per lattice 
site survive, and by identifying = S + = S x + iS y and n = S z + 1 the mapping to the spin 
| XXZ Hamiltonian is obtained. It reads 

h xxz = ~ £ 'S : r,,.s] -h^s;- t^sts* + sfs]} , (2) 

hO i (ij) 

where h = p — | ^oi- This Hamiltonian has been investigated extensively on the mean- 
field level for nearest-neighbor as well as next nearest-neighbor interactions. 

In the limit of a large number of bosons per lattice site one may parameterize = ^e l{f 
and take p equal to the boson density. Thus, only the phase remains as a dynamic variable 
and the Quantum-Phase Hamiltonian reads 

Hqp = \j2 n i u ij n j ~ A* n i ~ J cos&i ~ fj) , (3) 

hi i (ij) 

where J = pt. The phase and number are canonical conjugate variables and satisfy [tpi, rij] = 
iSij. This is exactly the Hamiltonian for Josephson junction arrays with an applied gate 
voltage in order to tune the chemical potential [ 2^| . It is convenient to introduce a parameter 
n o — ^/Y^iUij., which allows us to rewrite the Coulomb interaction and chemical potential 
term in Eqs. ([I]) and @ as 

H mt = \ J2( n i - no) u ij( n j ~ n o) ■ ( 4 ) 
A id 

We expect that the properties of the Quantum-Phase model will be periodic in the variable 
no with period 1. 
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III. MEAN-FIELD 



A first insight into the properties of these models may be gained by considering their 
mean-field solutions. The mean-field approximation for the spin model in Eq. (|2|) is straight- 



forward and was carried out by various authors [I I -13]. The question of a supersolid phase 
in this model was first addressed in Ref. [TI| and investigated further in Ref. ||13|| . We review 
the spin mean-field theory in Subsect. Ill A . A mean-field approximation for the Quantum- 
Phase model which is able to identify both diagonal LRO and off-diagonal LRO appeared 
in Ref. [H|. It is reviewed in Subsect. [IirB| Mean- field approaches to the Bose-Hubbard 
Model have been discussed in Refs. H0, but to our present knowledge there exists no 
mean-field approximation for the Bose-Hubbard model which allows for the identification of 
a supersolid phase. 

In Sect. 0we present a variational Ansatz for the many-body wave function that allows 
the three models to be treated on an equal footing. This is important, as it allows for 
a comparison between the models on the same level of approximation. Since we use a 
factorizable Ansatz for the wave function, it also has mean-field character. 

The validity of the mean-field approximation in two dimensions is often limited, since 
the lower critical dimension for models with continuous symmetries is two, which is the 
dimension we are particularly interested in. This is a serious problem at finite temperature 
where the models do not even exhibit superfluid long-range order as formulated in the 
theorem of Mermin and Wagner, and Hohenberg [|3T|. Useful results may, however, be gained 
from mean-field approximations at zero temperature where quantum dynamics formally adds 
an extra dimension, the imaginary time axis. The systems then do have long-range order 
and mean- field results may be qualitatively correct. This will be shown in the present article 
by scaling arguments and confirmed by the Monte Carlo simulations. In the remainder of 
this section we present the different mean-field type approaches. 



5 



A. Spin Mean-Field 



The spin mean-field theory can be formulated by linearizing the spin-spin interaction. 
This yields the standard self-consistency equations for the magnetization 



where the subscript i refers to the sublattice and (3 is the inverse temperature (we use 
units in which ks = h = c = 1). In order to identify spatially varying solutions one 
needs to introduce several sublattices, at least three for interactions including nearest and 
next-nearest neighbors. The symmetry in the xy-pl&ne reduces the number of independent 
variables and most of the phase boundaries can be determined analytically. If we allow for 
three different sublattices (A, B and C), the effective field is given by 



and corresponding equations for H B and He-. 

The connection to hard-core bosons is made by remembering that S z corresponds to the 
particle number, and S x ^ y to the boson creation operators. A staggered magnetization in in- 
direction implies a solid ( charge- density wave) ordering of the bosons. Finite magnetization 
in ^-direction translates into off-diagonal LRO for the boson model. The coexistence phase 
of both types of LRO is the supersolid phase for the bosons. The result is shown in Fig. [I] for 
U2 = 0.1 U\. Two different types of solid ordering occur for different chemical potential, i.e., 
(71", 7r) ordering around half-filling, and additional (n, 0) and (0, it) ordering around quarter 
filling. 




(5) 




(6) 



B. Quantum- Phase Model 



Now we turn to the mean-field treatment of soft-core bosons following Refs. flT4| , |32|| . 
After decoupling the hopping and the interactions between the different sites in Eq. (^), one 



arrives at a local problem 
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H* mf |^> = E i |^> , (7) 
with the 27r-periodic solutions for the wave function in the phase representation 

ipi{(Pi) = exp(ir]iXi) fi(xi) , (8) 

where (fi = 2 xt, r\i = 2 n — Uij{{nj) — n ). After these substitutions the following 

Mathieu equation for / is obtained 

r o 2 ^ 



dxf 



+ et + 2r i cos(2x i )\f i (x i ) = (9) 



4J 

with n = — ^2(cos((pi+p)) . (10) 

The sum runs over nearest neighbors f3 of i. The parameter is proportional to the energy 
and should be minimized. The ground state is defined via the self-consistency equation 

(nt) = I n t m 

(cos(^)) = I cos(^) |V°> ■ (11) 

In addition the free energy has to be minimized in order to have an unambiguous solution 
of Eq. (pTTj) . The numerical solution of this set of equations is straightforward. In order to 
allow for spatially varying solutions we again define sublattices. Our numerical solutions 



confirm the results of Ref . |L4| , however, some details differ. 

For on-site interaction we obtain Mott insulating lobes centered around integer values N 
of no, see Fig. |2|. In these incompressible lobes the density is pinned to N. At half integer 
values of no the critical value of the hopping vanishes, due to the degeneracy in energy of 
states with and A^ + 1 particles per site. 

Including nearest-neighbor interaction the calculation has to be performed on two sub- 
lattices. We obtain two different kinds of insulating lobes, one with integer filling and 
homogeneous density, another with half-integer filling and a checkerboard charge-density 
wave, or solid order, centered around half integer values of no- In addition we find a region 



of supersolid phase around the half lobe where the checkerboard charge- density wave and 
superfluidity coexist, see Fig. |3|. Again, the superfluid phase occupies the parts on the right 
of the phase diagram and also between the insulating lobes down to zero hopping, in slight 
contrast to the result of Ref. . 

The fraction of supersolid phase in the phase diagram depends on the ratio Ui/U , see 
Fig. |j. In the hard-core limit Uq — > oo the supersolid phase vanishes. An expansion of the 
ground state wave function |?/> t ) = J2n c n\ n ) m t° particle number eigenstates \n) yields more 
information about the nature of the supersolid phase. The inset of Fig. |] shows | C2 1 2 at 
no=0.5 as a function of J '/Uq. The figure indicates that double occupancy is important for 
the supersolid phase. Thus, we conclude that the supersolid phase is favored by fluctuations 
in the particle number. 

The inclusion of further interactions leads to more structure in the phase diagram. For 
next-nearest neighbor interaction we introduce three sublattices. As a result, The phase 
diagram shows a lobe with quarter filling and a phase where this quarter structure and 
superconductivity coexist, see Fig. |5|. 



C. Hard-core vs. Soft-core 



The previous subsections showed two extreme cases of the Bose-Hubbard model, the 
hard-core limit and the limit of large particle numbers. Although the phase diagrams look 
similar, these two limits behave qualitatively different as far as the supersolid phase is 
concerned. 

In the hard-core limit the existence of the supersolid is related to a finite next-nearest 
neighbor interaction. The supersolid does not exist in the model including nearest-neighbor 
interaction only ||13| . Furthermore, there is no supersolid phase at exactly half filling. This 
leads to the following interpretation of the supersolid phase. At densities corresponding to 
half filling the particles form an incompressible solid. Away from half filling superfluidity 
is enabled by defects in the solid structure that Bose-Einstein condense. This interpreta- 
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tion does, however, not explain why next-nearest neighbor interaction is necessary for the 
formation of the supersolid. 

In the limit of many particles per site the Quantum-Phase model applies. In this case 
the supersolid phase already exists for nearest-neighbor interaction only and at half filling. 
This is related to excitations which are forbidden in the hard-core limit. We observe that a 
large nearest-neighbor interaction U\ or small on-site interaction favors the supersolid, in 
the hard-core limit the supersolid is suppressed. In the soft-core limit the supersolid phase 
is also present at half-filling at the tip of the checkerboard lobe. Thus, it seems that the 
system itself generates the defects (particle-hole excitations) that Bose condense, thereby 
turning the solid into the supersolid. 



IV. VARIATIONAL METHOD 



Here we present a method which allows us to treat the three models within one 
scheme [27fl. The idea is to guess a variational wave function for the problem fi3f| . We 



use a variational wave function that is inspired by earlier successful Gutzwiller-like treat- 



ments of the spin model J34] and the Bose-Hubbard model p5| , p6[ . The occupation number 
representation is well suited for all three models. Explicitly we resolve the identity as 



i = EIW)(WI 



(12) 



and express the trace as 



Tr 



In,} 



E(WI-IW), 

{«;} 



(13) 



where the n.j are in [0, 1] for the spin | XXZ model, in [0..oo] for the Bose-Hubbard (BH) 
model, and in [—00. .00] for the Quantum-Phase (QP) model. Besides the difference in the 
allowed values of rii, the matrix elements differ for the three models, e.g., 



BH model: (Wi}\h\{ni}) = {{n' i¥lk }\{ni^ k ) (n' k \n k - l)v^fc 
QP model: ({n-}| exp itp k \{ni}) = (K^JIWfc) KK - l ) 



(14) 



As our variational wave function we take the product of single-site wave functions. This 
amounts to neglecting certain correlations and is therefore similar to a mean-field approxi- 
mation. We use the following Ansatz 



and kj and rrii are real variational parameters. The single-site wave function is a superpo- 
sition of boson number eigenstates weighted with a Gaussian. The average is controlled by 
rrii, the width by hi, and spatial variations (different sublattices i) are allowed. By mini- 
mizing the energy expectation value E = (ip\H\ip) with respect to k{ and rrii we obtain the 
ground state wave function within our approximation. In the following, expectation values 
are understood to be calculated with the ground state wave function. As the discrete sums 
can in general not be performed analytically, we treat the problem numerically. For on-site 
and nearest-neighbor interactions two sublattices arranged in a checkerboard configuration 
are introduced. We define the superconducting order parameter as (exp(z</?)) for the QP 
model, (b) for the BH model, and (S~) for the XXZ model. The structure factor 



yields information about the solid order. A non-vanishing structure factor at finite q sig- 
nals diagonal LRO. By decomposing the lattice into two sublattices we gain information 
about ,S(7r,7r). A finite S (71,71) corresponds to a checkerboard arrangement of the particles. 
For next-nearest neighbor interactions we introduce four sublattices which yields additional 
information about S(7t,0) and S(0,tt). 

Figure ^| shows the phase diagram for the QP model with on-site and nearest-neighbor 
interaction obtained with our variational Ansatz. Both the superfluid-insulator and the 
crystalline order transition can be identified and agree with previous mean-field calcula- 
tions |[4l ) p8[| . The phase boundaries are periodic in uq. We find that superfluidity sets in 





(15) 




(16) 
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simultaneously on both sublattices A and B. Within the supersolid phase the value of the 
superconducting order parameter on A differs from that on B, but both are nonzero. Beyond 
the crystalline order phase boundary the order parameter does not exhibit any difference on 
the sublattices. 

For the XXZ model the variational method confirms that there is no supersolid for on-site 
and nearest-neighbor interaction only. In the presence of next-nearest neighbor interaction 
we find supersolid phases in perfect agreement with |L3f and Fig. [l]. 



The phase diagram for the BH Hamiltonian is shown in Fig. [?]. The size of the lobes 
decreases with increasing n$. At point a the supersolid vanishes. This might be due to 
the lower bound for the occupation numbers n > 0. For small no charge fluctuations are 
suppressed. Hence, we conclude that charge fluctuations are necessary for the supersolid. 
For large no the phase boundaries of the QP model approach those of the BH model, if the 
hopping is rescaled according to Section [H], i.e., if t and J /no are identified. 

V. QUANTUM MONTE CARLO 

Several methods are available for determining the phase diagram for the three model 
Hamiltonians by numerical means. Exact diagonalization is possible for the XXZ model and 
small systems, i.e., less than about 26 lattice sites. Furthermore it is possible to perform 
Monte Carlo simulations for all the three models on much larger systems. The Bose- Hubbard 
model was studied in one dimension in Refs. []37| , j38 |, and one data point for the phase 



boundary was obtained in two dimensions in Ref. |p9||. The Bose-Hubbard model was more 



recently studied in the context of a supersolid phase in Ref. [20|. The two-dimensional 



Quantum-Phase model was studied in Ref. [(5J in relation to the universal conductivity at 



the SI transition and in Refs. [21, 28 1 in relation to supersolid phases 



In this section we present more results on Monte Carlo simulations of the two-dimensional 



QP model in addition to those in Ref. [28|. We discuss the on-site problem and compare to 



the analytic results oit/U expansions The inclusion of nearest-neighbor interactions 
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allows us to study the supersolid phase. 



A. Duality Transformation 

In order to study the Quantum-Phase model by means of Monte Carlo we map the 2D 
quantum model onto a (2+l)D classical model of discrete divergence-free currents. The 



essential feature of this mapping were presented in Ref. |23j and the derivation makes use of 



duality transformations that were developed in Ref. [43]. We start from the basic expression 



for the partition function Z for the Quantum-Phase model of Eq. (§) 

Z = Tr exp(-[3H QP ) , (17) 

where (3 is the inverse temperature. We go over to a Euclidean path-integral formulation by 
introducing time-slices, i.e., dividing f3 in iV T intervals of size e, such that N T e = j3. Inserting 
complete sets of states at each time slice we arrive at 

Z = / T^fhr ex P { - | Yl <'>: { ■,<'!■: + e f 1 Y *V 

{n iiT =0,±l,±2,---} ij,r i,r 

+ i Y n i^i,r + eJ Y cos^i.r - ¥j,r ~ Aj) } • (18) 

l,T (ij),T 

We included a coupling to a vector potential = (27r/$ ) // A-dl for later convenience. In 
order to proceed, we make one approximation, the Villain approximation for the cosine 
term. This amounts to expanding the cosine around all its minima 

exp{eJ cos(v3 ir - (p jT - A i:j )} w ex P{-— ^— ^(^r - W - 27rm ij>r - A^) 2 } , (19) 

where m is a directed discrete field that lives on the bonds between lattice sites. The function 
F is determined by demanding that the first two Fourier coefficients of the expressions on 
the left and right hand side of Eq. (19) are equal 0] and is given by 



F{X) 2x\og{I Q {x)/h{x)} • (20) 
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For large arguments F approaches 1, as is clear from a direct expansion of the cosine poten- 
tial. After a subsequent Poisson resummation, i.e., writing 

E fl m ioA = E J Vmf[m ijiT ]exp(2ni^2m ij:T J ijjT ) , 

{mij, T } {Jij.r} *?> 

an integration over the fields rriij jT and the phases </?j )T yields a representation in terms of 
divergence-free discrete current loops 

Z=[eJF(eJ)]- 1 ]T 5(V u J u )exp{-S[J}}, 
{^} 

j},T «,t \ / i,T,a=x,y i,r 

" "o) (*« + j^fo)) - «o) + E (^)i ( 21 ) 

ijY ir,a=x,y J i,r 

where v = x, y, r. The effective coupling constant is given by \J2/K with K = 8F(eJ)J/U , 
where e is defined implicitely by e = 1/ \JUoF(eJ)J. The time components J[ T of the current 
correspond to the boson numbers n^ T along their world lines. We keep track of the energy 
dependence of the determinant that arises from integrating out the m,ij jT . This is relevant 
for determining the energy density and specific heat by means of Monte Carlo. For zero 
magnetic field (A = 0) the action in the representation of Eq. Q2ID is real, and therefore 
suitable for Monte Carlo methods. Using the standard Metropolis algorithm we generate 
configurations of currents in a system of size L x L x L T . The condition that the currents 
J v be divergence-free is taken into account by making Monte Carlo steps that preserve the 
property V V J V = 0. Thus, we create or annihilate small current loops at every lattice site 
in all three directions, as well as periodic current loops that go through the whole system 
which is taken to have periodic boundary conditions in all three directions. The generation 
of configurations may be done in a canonical as well as in a grand-canonical way. Here 
we work in the grand-canonical ensemble in order to make contact to the mean-field phase 
diagrams. 

In the last line of Eq. (^T|) we used e = 1/ ^UoF(eJ)J. With this choice the couplings 
are isotropic for on-site interaction, leading to an efficient numerical algorithm. The choice 

13 




of the time lattice spacing e needs some justification. Introducing time slices in the partition 
function Eq. ([T^) is exact in the limit e —>■ 0. This would produce extremely anisotropic 
couplings in the current-loop model and the numerical simulation would become impossi- 
ble. Thus, a more reasonable choice for the time lattice spacing is necessary. The plasma 
frequency oo p = \/~JUo is a natural frequency for spin waves. A cutoff beyond this frequency 
by introducing a time spacing e ~ 1 /u p should not do any harm. We tested the dependence 
on the choice of e for the quantities of interest and observed that they are not sensitive to 
variations of e by one order of magnitude, provided that the temperature is kept constant, 
i.e., f3 = 1/T = const. = N T e. Additional justification for the choice of a finite time spacing 
e in the present paper is given by the fact that we are only interested in the critical regime, 
where a high-frequency cutoff should be irrelevant. 

B. Finite-Size Scaling 

The mapping to the current-loop model Eq. fl2~l"|) is well suited for simulating the behavior 
of the Quantum-Phase model with short range interactions. In practice we are able to study 
system sizes up to 12 x 12 x 12, or 10 x 10 x 25, respectively. The factor that limits the 
largest systems to have L= 12 is the necessity to make the periodic current loops in the 
spatial directions. Only they can change the superfluid stiffness, which is a topological 
quantity that cannot be changed by making local current loops. As phase transitions take 
place only in the thermodynamic limit, we need an additional ingredient that relates the 
data obtained on finite-size systems to infinite system size properties as the critical coupling 
and exponents. This is provided by finite-size scaling. 

Let us first consider the onset of superfluidity in the system. Therefore we study the 
behavior of the superfluid stiffness po, a quantity which measures the response to a twist in 
boundary conditions for the phase. A finite value of po in the thermodynamic limit reflects 
long-range phase coherence, i.e., off-diagonal LRO or superfluidity. A twist in the boundary 
conditions for the phases by an angle G increases the free energy density by 
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A/=^(!) 2 (22) 

The twist may be realized by applying a magnetic field. In general the frequency and 
wave-vector dependent stiffness p(k,ou IJi ) can be defined as 



,1\ 2 -5 2 lnZ 
pfcuj = — ^ 



2e7 ^ M ?^§ 



e-^T-tfcrj _ ^3) 

A=0 



The zero-frequency and zero wave- vector component of p(fc, u^) defines the superfluid stiff- 
ness. In terms of the currents it reads 

po = P(k=o^=o) = ^ (j2 J?)j "j = j^-H) ■ (24) 

The winding numbers w x = 4 ^ T ^ T can be measured easily; they are only modified by 
periodic current loops, local loops do not change their value. The symmetry in x/y-direction 
further simplifies the numerics. 

Direct measurement in the critical regime is ruled out by the divergence of the correlation 
length £ at the continuous transition. The powerful method of finite-size scaling, however, 
takes advantage of this fact and allows us to study the critical behavior. Near a continuous 
phase transition the diverging coherence length is cut off by the finite system size L. As 
a result all quantities will depend on the ratio £/L. A finite-size scaling relation for the 



superfluid stiffness po is readily derived [45| . Assuming hyperscaling, the critical part of the 
free energy density behaves as the inverse correlation volume £ d £ r - Here, £ is the correlation 
length in real space which may be different from the correlation length £ T in imaginary time, 
depending on the quantum dynamics of the model. The dynamical critical exponent z can 
be introduced through £ r oc £ z . Combining Eq. (p2|), the hyperscaling relation and the fact 
that £ is cut off by L (£ T by L T ) we arrive at the scaling behavior 

Po = L 2 ~ d - z p(bL l ^5,L T /L z ). (25) 

Here, b is a non-universal scale factor and p is a universal scaling function with a smooth 
dependence on L T /L Z and £/L = 5~ u / L. We assumed a power law critical behavior £ oc 5~~ u 
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with dimensionless distance to the transition 5 = (K — K*)/K* and critical coupling K*. At 
the critical point (5 = 0), L d+Z ~ 2 p is a function of L T /L Z only. Thus, plots of L d+Z ~ 2 p vs. 
K will intersect at the transition if L T /L Z is kept constant, with no further fitting involved. 
Furthermore, the data for L d+z ~ 2 po plotted as a function of L x l v b for different system sizes 
should collapse onto one single curve. This allows the exponent v to be determined. 

As we are interested also in diagonal LRO and a possible supersolid phase, we measure 
the structure factor as well. The structure factor is the (tt, 7r)-component of the equal time 
density-density correlation, see Eq. (0), and may be derived as the second variation of 
the free energy with respect to a staggered chemical potential. Expressed in terms of the 
currents J T it reads 

^ = TTr<E(-i) i+ ^;^>- (26) 

T ij,r 

The scaling relation for the structure factor is 

S n = L- wll, S{b'L 1/u 5, L T /L Z ) . (27) 

This scaling relation arises, since the structure factor S n is related to the square of the 
staggered magnetization order parameter M = Yl,i{~^-) 1 <J[ for which the exponent (3 is the 
order parameter exponent. Near the phase transition, M ~ 5@ and £ ~ 5~ u , from which the 
quoted form follows. From a three-parameter fit to the scaling relation for different system 
sizes, the exponents v and (3 as well as the critical coupling constant K* are determined. 
Again, plots of L 2l3 / u S n will cross at the critical point if the ratio L T /L Z is kept constant. 

In both scaling relations a knowledge of z is assumed. This is well-known for the SI 
transition and the scaling of po, i.e., in general z—2, except for the tips of the lobes where 
z—1 as dictated by particle hole-symmetry. This argument is based on a coarse-graining 
treatment || that allows for the explicit construction of a Ginzburg-Landau- Wilson free 



energy functional for the superconducting order parameter P,^6|-|4q|, from which z can be 
read off. For the transition related to the structure factor the situation is less clear. We were, 
for instance, not yet successful in explicitly constructing the Ginzburg-Landau- Wilson free 
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energy functional for the CDW transition. From symmetry arguments one expects z—1 at 
the mean- field level |fI5 |. Also the spin- wave analysis of Ref. |2(J predicts z—1. However, it 
is unclear whether or not z for the CDW will be renormalized by the coupling to the gapless 
sound mode in the superfluid part. As a working hypothesis we take z=l and comment later 
on the consistency of this assumption. 

C. Numerical Results 

We report about the SI transitions for on-site interaction and nearest-neighbor interaction 
in the complete range of the chemical potential. A phase diagram is mapped out and the 
scaling predictions are verified. The existence of a supersolid phase in a finite region of the 
phase diagram is firmly established. 

The case of on-site interaction has been extensively studied in the past years. A phase 
diagram consisting of lobe-like structures in the t — p or J — n plane is found. Mott- 
insulating lobes are centered around integer values of Uq. This has been observed in mean- 



field approximations, t/U expansions and in Monte Carlo simulations of the one- 



dimensional Bose Hubbard model |37| . The shape of the lobes in two dimensions is yet 
unknown. Therefore we first map out this phase diagram which is shown in Fig. |^. We 
performed the simulations for fixed n Q and tuned through the transition by varying the 
effective coupling K. In the phase diagram in Fig. ||] this corresponds to moving on horizontal 
lines through the phase transition. 

At integer values of uq the system is particle-hole symmetric and we have z = 1 . Accord- 
ing to the scaling form of Eq. ([25|) we use L = L T which keeps the second argument of the 
scaling function constant. Four different system sizes with L — 6,8, 10, 12 where studied. 
For the largest systems we typically use 10 6 sweeps through the system for equilibration 
and measurement. In each sweep we try all possibilities of local and periodic current loops. 
Without any fitting the curves of poL d+z ~ 2 = poL vs. K should cross in the critical point K* . 
This serves as a good test of the scaling prediction. The critical exponent v can be extracted 
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by fitting the data of different system sizes near the transition to one scaling curve. This is 
shown in Fig. [8] where we plot the data as a function of 5L l / u . For n = we obtain the 
critical coupling K* = 0.886 ± 0.003, and using a x 2 fit ^ = 0.69 ± 0.06. The exponent v 
agrees with the expected 3D-XY behavior, where v « 2/3. 

Away from the symmetry points at the tips of the lobes the dynamical critical exponent 
is z = 2. In order to keep L T /L Z constant we simulate systems with L T = L 2 /4, i.e., 
LxLxL T = 6x6x9, 8x8x16 and 10 x 10 x 25. Now the curves of p L 2 should cross in 
the critical point. We obtained data for hq =0.05, 0.1, 0.45 as shown in Table |. From a 
fit to the scaling form we find 0.42 < v < 0.52 with errors of ~ 20%. This is in agreement 
with the expected mean- field transition where v = 1/2. A more accurate determination of 
the critical exponents turns out to require an enormous numerical effort with our method 
and is beyond the scope of our present work. 

In Fig. 0, the Monte Carlo data are shown together with mean-field phase boundaries 
and the results of the t/[/-expansion of Refs. |4T|J42]] in the limit of large n Q (or //). A lobe 
shape is observed, these lobes are sharper than predicted by mean-field, but smoother than 
expected from the t/f/-expansion. Approaching no = 0.5 we observe that the critical value 
of J/Uq decreases and a simple extrapolation yields J* = for n = 0.5. A numerical answer 
to the question whether J* is zero at no = 0.5 is not possible, since the acceptance rates are 
very small for small J/Uq. 

The inclusion of nearest-neighbor interactions yields a richer phase diagram. The scaling 
for the superfluid stiffness is not modified by the inclusion of short-range interactions, i.e., the 
universality class is preserved. We expect z = 1 in the case of particle-hole symmetry, z = 2 
otherwise. For nearest-neighbor interaction the lines with integer and half-integer values of 
n exhibit particle-hole symmetry. For the transition related to charge- density wave order, 
we expect z — 1 throughout the phase diagram. Guided by the mean-field phase diagram 
we expect the existence of two different kinds of lobes, with homogeneous particle densities 
centered around integer values of n , and with a checkerboard charge-density wave centered 
around half integer values of n . The supersolid phase is expected in the vicinity of the 
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checkerboard lobe. 

For n = and 0.5 we simulated Lx Lx L systems, where L— 6, 8, 10, 12, as suggested by 
particle-hole symmetry and z — 1. For rio=0.1, 0.2 and 0.4 we obtain z = 2 for the superffuid 
transition. In order to keep the ratio L T /L Z constant, we simulated L x L x L 2 /4 systems, 
where L= 6, 8, 10. For n =0.4 near the transition for the diagonal LRO we performed the 
calculations with both z = 1 and z — 2. The results are summarized in Fig. |9| and Fig. 
and Table [n[ In the following we discuss our data in detail. 

First we present our data for no=0.5. Figure |9] shows that there are two separate tran- 
sitions for diagonal and off-diagonal LRO with a coexistence region in between where both 
the superffuid density and the structure factor scale to a finite value in the thermodynamic 
limit. This demonstrates the coexistence of diagonal LRO and off diagonal LRO for soft-core 
bosons with nearest-neighbor interaction in 2 dimensions. Table shows that the exponent 
v is different for the two transitions. For the transition related to superfluidity (point j3 in 
Fig. H|) we find v = 0.65 ± 0.08 which is consistent with the 3D XY universality class. For 
the transition related to crystalline order (point 7) the universality class is not known. We 
find v ^0.55 and j3 ^0.21. 

Also at ?2o=0.4 we find two separate transitions that are the boundaries for the supersolid 
phase in between, see Table |T[ As compared to uq=0.5 both transitions are shifted to 
smaller values of the coupling constant K. This is consistent with the mean-field phase 
diagram. Again the two transitions have different critical exponents. The transition related 
to superfluidity (the line separating phases "Sol" and "SSol" in Fig. |3|) has v = 0.44 ± 0.08 
which is consistent with a mean-field transition in d + z = 4 effective dimensions. For the 
transition related to crystalline order (between phases "SSol" and "SF") we compare the 
data for z = 2 and z — I 'm Fig. [10]. A fit to the expected scaling behavior is equally possible 
in both cases. The statistical errors of the data for z = 2 are larger, since the simulation 
for z = 2 requires larger L T (up to 25) which decreases the acceptance ratio for periodic 
loops in this direction. The values for the critical coupling and the critical exponents for 
both z = 1 and z = 2 are shown in Table |J, they coincide within the error bars. With the 
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values of Table [IT] the scaling relation 2(3 = v[d + z — 2 + 77) leads to 0.8 ± 0.2 = z + rj, with 
correlation function exponent 77 (small and usually positive). This rules out z = 2 and one 
is led to conclude that indeed z ~ 1 for the crystalline transition, independent of no- 

The small value of (3 rules out mean-field behavior for the CDW transition. Insight into 
the nature of this transition is gained by the following consideration. In the neighborhood 
of this transition, i.e., far in the superfluid phase, the x,y-components of the currents J 
fluctuate strongly and may be integrated out as Gaussian fluctuations. In other words: one 
expects the staggered magnetization order parameter to couple to the gapless 4 th sound 
mode of the superfluid. This yields strong long-range interactions in the time direction for 
the r-components of the currents J . It is likely that these long-range interactions are a 
relevant perturbation and suppress the exponent (3. 

Finally, we discuss the data for no = 0, 0.1, and 0.2. In these cases there is only one phase 
transition, as the Mott-insulating lobes (phase "MI" in Fig. |3|) do not have any non-trivial 
crystalline order. Our data are consistent with a transition in the 3D XY universality class 
for Uq=0 and a mean-field phase transition for no = 0.1 and 0.2. 

VI. DISCUSSION 

We studied the T = phase diagram of two-dimensional interacting bosons by various 
techniques. The combination of the various aspects of these treatments allow us to draw the 
following conclusions. 

For on-site interaction we have seen the lobe structure of the phase boundary in the 
t — fi or J — n plane. This structure is exactly periodic in the chemical potential no in the 
Quantum-Phase model of Eq. (^|). In mean field or variational treatments of the different 
models, the lobes are parabolic near their tips. The Monte Carlo results, however, show 
that the lobes are sharper than predicted by mean-field approximations. This deviation 
near the tips of the lobes is expected by scaling arguments: due to the additional particle- 
hole symmetry at the tips of the lobes, the effective dimension is reduced to d e ff = d + z = 3 
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and enhanced fluctuations destroy part of the ordered (superfluid) phase. This is consistent 
with the value of the critical exponent v obtained from the simulations [y sa 2/3) that agrees 
well with a 3D XY critical point. Similar sharpening of the lobes is seen in a higher-order 



t/U expansion for the BH model [O 



A central point of this paper is the consideration of finite-range interactions. They lead 
to a richer phase diagram, due to commensurability effects. In particular, charge-density 
wave or solid order appears in parts of the phase diagram. In our definition the solid phase 
is characterized by the breaking of the translation symmetry of the underlying lattice, i.e., 
if a density wave with a wave vector smaller than the lattice is present. 

Upon inclusion of nearest-neighbor interactions new half-filled lobes around half-integer 
n appear in the phase diagram, in which the checkerboard (tt, it) density wave is stable. 
For next-nearest neighbor interactions additional density waves exist in which 1 or 3 out of 
4 sites are occupied, around no— 1/4 or 3/4. These pure density waves are incompressible 
Mott-insulators. Surprisingly, each Mott-insulating charge-density wave seems to have a 
corresponding compressible supersolid phase in which the particular CDW order coexists 
with off diagonal LRO. 

For soft-core bosons the supersolid phases appear already for nearest-neighbor inter- 
action. Hard-core bosons, on the other hand, supersolidify only if at least next-nearest 
neighbor interaction is present. Another difference is the presence of a supersolid at the 
tip of the half lobe: hard-core bosons only form a supersolid away from half-integer filling, 
whereas soft-core bosons do so also at half-integer filling. 

The difference between Bose-Hubbard and Quantum-Phase models is another issue. Our 
results indicate that, on the mean-field level, the two models are, apart from a trivial rescal- 
ing of the hopping matrix element, almost identical for large fi or n . For small fi corre- 
sponding to a filling of less than one boson per site the mean-field solution for the two models 
differ essentially and are not described by a simple rescaling of the parameters. Within our 
mean-field approximation we find that the supersolid phase is suppressed around the solid 
lobe with exactly half filling (alternating empty and singly-occupied sites). This may be one 
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of the reasons why in a recent numerical simulation of the BH model [20| a supersolid phase 
was not found for exactly half filling. 

Our combined mean-field and Monte Carlo analysis has shown that the density-wave 
order transition is not mean-field like. Both the exponents (j3 ~ 0.21 and v « 0.55) and the 
location of the phase boundary deviate considerably from the mean-field results (/3 = 1/2 
and v = 1/2), see also Fig. [| This may be related to the coupling of the CDW order 
parameter to the gapless 4 th sound mode in the superfluid and remains a subject for further 
study. A related issue is the value of the dynamical critical exponent z for this transition. 
The value z—1 used in the simulations gives good scaling. However, other values of z turn 
out to give reasonable scaling as well and yield identical values for /3, v and K*. It seems that 
the CDW transition is rather indifferent to the value of z used to scale the systems studied 
numerically. From the scaling relation 2/3 jv = d + z — 2 + 77 we expect the combination 
z+i] « 0.8±0.2 which would indicate either a small rj and renormalized z ~ 0.8 or z ~ 1 and 
r) negative. Further simulations and renormalization group studies are required for solving 
this problem. 

In summary we performed an analysis of interacting bosons in 2D within zero tempera- 
ture lattice models, using and comparing both mean-field theory and exact Quantum Monte 
Carlo simulations. We obtained the full phase diagram for on-site interaction as well as 
for nearest-neighbor interactions, in which case our simulations establish the existence of a 
supersolid phase in which a CDW and superfluidity coexist. 

Acknowledgments: We thank G. T. Zimanyi for motivation and discussions, and J. K. 
Freericks and H. Monien for discussions and their data as plotted in Fig. |2j. The work 
was partially supported by 'Sonderforschungsbereich 195' of the DFG and by the Swiss 
Nationalfonds (AvO). 
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TABLES 



n 


z 




J*/U 


V 


n 


z 


K* 


r/u 


V 


0.00 


1 


0.886 ±0.003 


0.152 


0.69 ±0.06 


0.25 


2 


0.64 ±0.01 


0.098 


0.5 ±0.1 


0.05 


2 


0.85 ±0.01 


0.144 


0.5 ±0.1 


0.30 


2 


0.56 ±0.01 


0.081 


u 


0.10 


2 


0.82 ±0.02 


0.137 


u 


0.35 


2 


0.47 ±0.01 


0.062 


u 


0.15 


2 


0.77 ±0.02 


0.126 


u 


0.40 


2 


0.37 ±0.01 


0.042 


u 


0.20 


2 


0.71 ±0.01 


0.113 


u 


0.45 


2 


0.26 ±0.01 


0.023 


u 



TABLE I. Critical couplings and exponents for on-site interactions 



the transition for 




off-diagonal LRO 


diagonal LRO 


n 


z 


K* 


J*/U 


V 


z 


K* 


J*/U 


V 




0.5 


1 


0.775 ± 0.005 


0.127 


0.65 ±0.08 


1 


0.837 ±0.005 


0.141 


0.55 ±0.05 


0.21 ±0.04 


0.4 


2 


0.645 ± 0.008 


0.099 


0.44 ±0.08 


2 


0.749 ± 0.006 


0.122 


0.5 ±0.1 


0.25 ±0.10 


0.4 










1 


0.747 ±0.007 


0.121 


0.55 ±0.06 


0.21 ±0.08 


0.2 


2 


0.446 ± 0.005 


0.057 


0.5 ±0.1 


for comparison: 


0.1 


2 


0.707 ±0.007 


0.112 


0.49 ±0.11 


mean-field: 


1/2 


1/2 


0.0 


1 


0.843 ± 0.005 


0.142 


0.61 ±0.08 


3D XY: 


2/3 


1/3 



TABLE II. Critical couplings and exponents for the different transitions with nearest-neighbor 
interactions U\/Uq = 0.2 
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FIGURES 

FIG. 1. : Phase diagram for hard-core bosons, as obtained from spin mean-field theory on 
the XXZ model. Nearest and next-nearest neighbor (U2/U1 = 0.1) interactions are included. 
Several phases appear: superfluid (SF), Mott-insulating (MI), checkerboard solid (Sol) and the 
corresponding supersolid (SSoll), and quarter-filled (QF) solid with the corresponding supersolid 
(SSol2). 

FIG. 2. : Phase diagrams for soft-core bosons, as obtained from the mean-field analysis of the 
QP model with on-site interaction Uq only (round curves on the left). The symbols are the Monte 
Carlo data as discussed in Sect. The sharp lobes are results of the third-order i/[/-expansion 
of Ref. [H^] (extreme right) and the extrapolation to infinite order in t/U from Ref. |42j. The 
superfluid phase is denoted by "SF" and the Mott-insulating phase by "MI" . 

FIG. 3. : Phase diagrams for soft-core bosons, as obtained from the mean-field analysis of 
the QP model with on-site (Uq) and nearest-neighbor (U\/Uq = 1/5) interaction. The symbols 
are the Monte Carlo data as discussed in Sect. [v|. The checkerboard charge-density wave is de- 
noted by "Sol", the supersolid phase by "SSol", the superfluid phase is denoted by "SF" and the 
Mott-insulating phase by "MI" . 

FIG. 4. : Supersolid region "SSol" at uq= 0.5 as a function of U\/Uq in the mean-field ap- 
proximation of Section [II E| . Inset: Occupation-number probability | C2 1 2 at uq= 0.5 for the two 



sublattices A and B at the particular value of U\/Uq = 0.2. 

FIG. 5. : Mean-field phase diagram for nearest-neighbor and next-nearest neighbor interac- 
tions, Ui/Uq = 0.2, U2/U0 = 0.02. In addition to the phase diagram in Fig. ||, small lobes at 
quarter filling and tree-quarter filling appear at small J with their corresponding supersolids. The 
phases are labeled as in Fig. [j] . 

FIG. 6. : Phase diagram for the Quantum-Phase model as determined from the variational 
Ansatz with nearest-neighbor interactions, Ui/Uq = 1/5. 
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FIG. 7. : Phase diagram for the Bose-Hubbard model as determined from the variational 
Ansatz with nearest-neighbor interactions, U\/Uq = 1/5. For small uq the supersolid phase vanishes 
at point a. 

FIG. 8. : Scaling plot of the Monte Carlo data for the superfluid stiffness po for no = with 
on-site interaction only. 

FIG. 9. : Monte Carlo data for the superfluid stiffness po and the structure factor S n for 
no = 0.5 with nearest-neighbor interaction U\/Uq = 1/5. We clearly observe two distinct transitions 
at the crossing points of these scaling plots, and a coexistence phase in between. 

FIG. 10. : Scaling plots for the structure factor S w for z = 1 (lower curve) and z = 2 (upper 
curve) at no= 0.4 and U\/Uq = 1/5. 
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